
__author__="martin"
__date__ ="$18.12.2009 19:05:18$"

import scipy
import scipy.ndimage
import pylab
from PyQt4 import QtGui, QtCore
import sys

def filter(mat, f, parent=None):
    xsize, ysize= mat.shape
    xrange=range(1, xsize-1)
    yrange=range(1, ysize-1)
    max_n=xsize*ysize
    mat2=scipy.zeros([xsize, ysize])
    mat3=scipy.zeros([3, 3])
    n=0
    if parent != None:
        parent.generate_progress_dialog("filtering...", "Cancel", max_n)
    for x in xrange:
        for y in yrange:
            n+=1
            if parent != None:
                parent.show_progress(n)
            mat3=mat[x-1:x+2, y-1:y+2]
            mat3=mat3*f
            mat2[x ,y]=sum(sum(mat3))
    if parent != None:
        parent.close_progress_dialog()
    return mat2

def create_matrix_from_point_list(point_list, xwidth, ywidth, parent=None):
    mat=scipy.zeros([ywidth, xwidth])
    p_array=scipy.array(point_list)
    #print p_array[:, 0]
    x_min=p_array[:,0].min()
    x_max=p_array[:,0].max()
    y_min=p_array[:,1].min()
    y_max=p_array[:,1].max()

    dx=abs(x_max-x_min)
    if dx==0.0:
        dx=1.0
    dy=abs(y_max-y_min)
    if dy==0.0:
        dy=1.0
    print dx, x_min, x_max
    p_range=range(0, len(p_array))
    if parent != None:
        parent.generate_progress_dialog("working...", "Cancel", len(p_array))
    #progress=QtGui.QProgressDialog("Working...", "Cancel", 0, len(p_array) )
    for n in p_range:
        #progress.setValue(n)
        if parent != None:
            parent.show_progress(n)

        x=int((p_array[n,0]-x_min)/dx*xwidth)
        y=int((p_array[n,1]-y_min)/dy*ywidth)
        #print x, y
        if 0<=x<xwidth and 0<=y<ywidth:
            mat[y, x]+=1
    if parent != None:
        parent.close_progress_dialog()
    return mat

def make_ifs_pic(xwidth=800, ywidth=600, max_iterations=10**6, function_system_filename="", transformation_functions_filename="", gamma=1.2, parent=None):

    def iterate( x, y, f):


        j=int(scipy.random.randint(len(f)))
        #print x, y
        #print j
        try:
            (nx, ny)=eval(f[j])
        except:
            (nx, ny)=(x, y)
        #print x, y
        return (nx, ny, j)

    def load_transformation_functions(transformation_functions_filename):
        f=list()
        try:
            file=open(transformation_functions_filename)
        except:
            print "file not found!"
            return f
        for line in file:
            f.append(line.strip())
        file.close()
        return f

    def transform(x, y, function_list):
        phi=scipy.arctan(x/(y+10**(-6)))
        omega=scipy.arctan(y/(x+10**(-6)))
        #print x, y
        try:
            r=(x**2+y**2)**0.5
        except:
            r=1
        psi_1=scipy.random.random()
        psi_2=scipy.random.random()
        #print w
        x_t=0.0
        y_t=0.0
        x_f=0
        y_f=0
        #print function_list
        for f in function_list:
            try:
                x_f, y_f=eval(f)
            except:
                x_f, y_f= x, y
            if (scipy.isfinite(x_f)):
                x_t+=x_f
            if (scipy.isfinite(y_f)):
                y_t+=y_f
        return x_t, y_t

    def load_function_system(function_system_filename):
        f=list([])
        try:
            file=open(function_system_filename)
        except:
            print "file not found!"
            return f
        for line in file:
            data=line.strip().split(";")
            n_range=range(0, int(float(data[1])*100))
            for n in n_range:
                f.append(data[0])
        file.close()
        return f


    x, y= 1, 1
    #max_iterations=10**6
    i=0
    old_pc=0
    p_list=list([])
    function_system=load_function_system(function_system_filename)
    transform_function=load_transformation_functions(transformation_functions_filename)
    if parent != None:
        parent.generate_progress_dialog("working...", "Cancel", max_iterations)
    # progress=QtGui.QProgressDialog("Working...", "Cancel", 0, max_iterations )
    print function_system
    while i <= max_iterations:
            i+=1
            if parent != None:
                parent.show_progress(i)
            #progress.setValue(i)
            x, y, c=iterate(x, y, function_system)
            x, y=transform(x, y, transform_function)
            xt, yt=x, y
            p_list.append([xt, yt, c])
  
    if parent != None:
        parent.close_progress_dialog()
    mat=create_matrix_from_point_list(p_list, xwidth, ywidth)

    mat=scipy.log(mat+1)/scipy.log(mat.max()+1)
    mat=mat**gamma
    #print mat
    #mat=mat/float(mat.max())
    f=scipy.zeros([3,3])
    f[0:3, 0:3]=1/16.0
    f[1, 1]=8/16.0
    #mat=filter(mat, f, parent)
    mat=scipy.ndimage.convolve(mat, f)


    #mat=mat/float(mat.max())
    return mat


def make_julia_pic(xwidth=800, ywidth=600, c=-1, max_iterations=10, parent=None):

    def iterate(number_of_iterations, value, x0, y0, c=0):
        z=complex(x0, y0)
        n=0

        while (n < number_of_iterations) and (abs(z) <= value):
            z=z**2+c
            n+=1
        return n

    mat=scipy.zeros([ywidth, xwidth])
    xrange=range(0, xwidth)
    yrange=range(0, ywidth)
    max_points=xwidth*ywidth
    i=0
    if parent != None:
        parent.generate_progress_dialog("working...", "Cancel", max_points)
    for x in xrange:
        #print x
        for y in yrange:
            i+=1
            if parent != None:
                parent.show_progress(i)
            mat[y, x]=iterate(max_iterations, 2, (x-xwidth/2)/float(xwidth)*4, (y-ywidth/2)/float(ywidth)*4, c)

    if parent != None:
        parent.close_progress_dialog()
    return mat

def make_curlicue_pic(xwidth=800, ywidth=600, s=(1+5**0.5)/2.0, gamma=1.2, max_iterations=10**6, parent=None):
    mat=scipy.zeros([ywidth, xwidth])
    i=0
    omega=0
    phi=0
    x, y=0, 0
    old_pc=0
    e=scipy.e
    pi=scipy.pi
    p_list=list()
    max_value=0
    if parent != None:
        parent.generate_progress_dialog("working...", "Cancel", max_iterations)
    while i <= max_iterations:
        i+=1
        if i*100/max_iterations > old_pc:
            old_pc=i*100/max_iterations
            if parent != None:
                parent.show_progress(i)

        omega=(omega+2*pi*s)%(2*pi)
        phi=(phi+omega)%(2*pi)
        x, y=x+scipy.cos(phi), y+scipy.sin(phi)
        p_list.append([x, y,0])

    mat=create_matrix_from_point_list(p_list, xwidth, ywidth)

    mat=scipy.log(1+mat)/scipy.log(2)
    mat=mat**gamma
    if parent != None:
        parent.close_progress_dialog()
    return mat

def make_newton_pic(xwidth=800, ywidth=600, n=3):

    def iterate(z0, n):
        z=z0
        z_old=0
        #n=0
        epsilon=10**(-3)
        #print z, n
        while (abs(z-z_old) > epsilon ) :
            z_old=z
            if z!=0:
                f_z=(n*z**(n-1))
            else:
                f_z=0
            if f_z != 0:
                z=z-((z**n)-1)/f_z
            else:
                z=0

        return z

    mat=scipy.zeros([ywidth, xwidth])
    xrange=range(0, xwidth)
    yrange=range(0, ywidth)

    for x in xrange:
        print x
        for y in yrange:
            z=complex((x-xwidth/2.0)/float(xwidth/2.0)*3, (y-ywidth/2.0)/float(ywidth/2.0)*3)
            i_z=iterate(z, n)
            mat[y, x]=i_z.real*i_z.imag

    mat=mat-mat.min()
    return mat

if __name__ == "__main__":
    crange=pylab.linspace(-1, 1, 100)
    app=PyQt4.QtGui.QApplication(sys.argv)



    mat=make_ifs_pic(400, 400, "farn.ifs")
    #mat=mat-mat.min()
    #mat=scipy.log(1+mat)/scipy.log(2)
    #mat=mat/float(mat.max())
    #print mat.mean()
    pylab.matshow(mat)
    pylab.gray()
    pylab.show()
    sys.exit(app.exec_())

    
